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Oh ; Abstract 

Propagation of extremely short electromagnetic pulses in a homogeneous doubly- 
resonant medium is considered in the framework of the total Maxwell-Duffing-Lorentz 
model, where the Duffing oscillators (anharmonic oscillators with cubic nonlineari- 
ties) represent the dielectric response of the medium, and the Lorentz harmonic os- 
■ cillators represent the magnetic response. The wave propagation is governed by the 

one-dimensional Maxwell equations. 

It is shown that the model possesses a one-parameter family of traveling-wave solu- 
"sj- ■ tions with the structure of single or multiple humps. Solutions are parametrized by the 

£S) ■ velocity of propagation. The spectrum of possible velocities is shown to be continuous 

on a small interval at the lower end of the spectrum; elsewhere the velocities form a dis- 
crete set. A correlation between the number of humps and the velocity is established. 
The traveling-wave solutions are found to be stable with respect to weak perturbations. 
Numerical simulations demonstrate that the traveling-wave pulses collide in a nearly 
^ ' elastic fashion. 
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1 Introduction 

The recent demonstration of artificial materials (metamaterials) with the left oriented triplet 
of electric E, magnetic H and wave vector k of electromagnetic field El, HI, 3, H| stimulated 
studies of nonlinear optical phenomena in such materials @, 0, @, 0, LLO, U|. Nonlinear 
dynamics of extremely short optical pulses in left-handed materials was the subject of partic- 
ular interest in several recent papers 0, 0, 13, 14]. The first experimental realization of the 



left-handed property based on the resonant response of the artificial material to both elec- 
tric and magnetic fields was described in [l|. To mention just one of the latest experimental 
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achievements, Valentine et al 0] were able to observe the negative refractive index in the 
balk material in the optical range. A theoretical description of the electromagnetic wave in- 



teraction with such double resonance materials (DRM) was considered in [15|, [ly, |17J, UM, [19 



Presence of two frequency intervals with different orientation of (E, H, k) triplets is a char- 
acteristic feature of such materials. 

Most of the studies of electromagnetic pulse propagation in DRM has been conducted 
in the slowly varying envelope approximation. On the other hand, there is a broad area 
of nonlinear optical phenomena taking place in the limit of extremely short pulses, when 



the slowly varying envelope approximation is not valid 20]. The case of extremely short 
electromagnetic pulses offers a new type of nonlinear interaction, when different frequency 
components of electromagnetic pulses have different orientations of the (E, H, k) triplets. 

The design of currently available DRM is based upon the use of embedded metallic 
structures whose size is on the same order as the spatial size of an extremely short electro- 
magnetic pulse. Therefore a theoretical and numerical investigation of the currently existing 
DRM would require 3D computer simulation on Maxwell's equations that takes into account 
the strong inhomogeneity of composite materials. Recently, there have been introduced 
some qualitatively different approaches to design of DRM, including the use of multilevel 
atoms 2l|, |24j, |25|; the latter gives rise to a spatially homogeneous medium. Possibilities 
of experimental realizations of such an approach were recently discussed in 22|, |23[ . As a 



first step in the theoretical investigation of electrodynamics of homogeneous DRM in this 
paper we study a simple model of a homogeneous doubly-resonant medium. Even under 
such simplification, dynamics of extremely short pulses turn out to be quite complex. 



2 Basic equations 

The system of equations that describe interaction of coherent light with a medium consisting 
of molecules (considered as harmonic oscillators) is known as the Maxwell-Lorentz model [271 ] . 
In this work we use a version of the Maxwell-Lorentz system that is extended to account 
for simultaneous magnetic and electric resonances, with the magnetic susceptibility being 
linear, while the electric polarization being nonlinear. Consider the general form Maxwell's 
equations: 

V x E = -c' x B u V x H = -c~ l D t (2.1) 
B = H + 4vrM, D = E + 4ttP 

For simplicity, we consider transverse electromagnetic plane waves propagating along the 
z-axis with the electric field E = (E(z, t), 0, 0) and the magnetic field B = (0, B(z, t), 0). 
Then the Maxwell equations transform to the scalar form: 

dE + 1 dB _ Q dH + 1 dD _ Q 

dz c dt ' dz c dt 

B = H + 4nM, D = E + A-kP (2.3) 
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which leads to 

E z + c~ l H t = -Anc-Hlt, H z + c^Et = -^KC~ x P t (2.4) 

The system (12 .4p must be closed by two additional equations describing the interaction of 
the electric and magnetic fields with the DR medium. As usual, it is convenient to avoid the 
47r-f actors by changing the units for M and P: M = 4irM, P = Air P. In the sequel we drop 
the tildes over M and P. 

Assume that the medium polarization is defined by the plasma oscillation of electron 
density, 

Ptt = ^ 2 p E 

Here u p is an effective parameter characterizing polarizability of the medium; in the case 
of metallic nanostructures it would be the effective plasma frequency. To account for the 
dimensional quantization due to the confinement of the plasma in nanostructures one should 
include the additional term LVpP, where uj£, is the frequency of dimensional quantization. 
We take into account nonlinearity in the lowest order of P, which is P 3 . A more accurate 
analysis, based on a quantum mechanical approach 29] and experimental measurements 3ol 



confirms validity of this assumption. Therefore we consider the modeling equation for the 
medium polarization dynamics in the following form 

P tt + cu 2 D P + nP 3 = cu 2 p E (2.5) 

where k is a constant of anharmonisity. To account for magnetic resonances we use the 
standard model [151 

M tt + u 2 T M = -(3H tt (2.6) 

Here j3 is a constant characterizing magnetization. 

We represent equations (12.41) . (12.51) and (12.61) in a dimensionless form by introducing 
t = t/r (t = 1/uj p is the characteristic time), r\ = z/z (z = cr is the characteristic 
distance), q = P/Pq (Pq = oj p /^fhl is the maximal achievable medium polarization). It is 
convenient to normalize remaining variables as follows: m = M/Pq, e = E/Pq, h = H/Pq. 
The system of dimensionless equations then takes the following form: 

h/f I — iTi^-^ 

q TT + uj\q + 7 g 3 = e (2.7) 
m TT + oj 2 m = —j3h TT , 

where 7 = k/ lo x = uj d /uj p , u> 2 = ujt/oj p . 

The system possesses the following conserved quantity: 

~ J [fiul (e 2 + uj\q 2 + |g 4 ) + (3uo 2 2 (h + mf + u 2 (1 - (3) m 2 (2.8) 
+(3uj 2 2 (q r f + [m T + f3h T ] 2 ] dr] = 
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which is positive-definite for (3 < 1. For the traveling- wave solutions the conservation relation 
(I2.8P yields conservation of electromagnetic energy 

-J (e 2 + /i 2 ) drj = const 

(see [28] for details). A natural question arises is whether the system in (12.71) possesses any 
solitary-wave solutions. We address this issue in the following section. 



3 Solitary wave solutions 

Consider a traveling wave solution of (12. 7p . i.e., a solution that is a function of the variable 
( = t — r)/V. Then the PDEs in (12.71) become ODEs, and one obtains the following system: 



h' - e'/V = -m (3.9) 

e' - h'/V = -q' (3.10) 

q" + ulq + 7 g 3 = e (3.11) 

m" + u 2 2 m = -/3h" (3.12) 



Upon the integration of equations (13.91) and (13.101) once, we get the algebraic conservation 
relations 

Vh-e = -mV + R 
-h + eV = -qV + S 

We are interested in a traveling-wave solution on the zero background, hence h = m = q = 
e = at ±oo; therefore the constants of integration R — S — 0. This yields the following 
expressions for h and e 

h = aim + a 2 q (3.13) 
e = a<im + a\q (3-14) 

where 

a x = V 2 (1 - V 2 ) -1 , a 2 = V (1 - V 2 )' 1 (3.15) 

We insert expressions (13.131) and (13. 14f) for h and e into the equations (13. lip and (13.121) 
for q and m and obtain the following system of second order equations: 

q" + (ujI - a x ) q - a 2 m + 7g 3 = 
(3a 2 q" + (1 + j3a x ) m" + u 2 2 m = 
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This system can be diagonalized with respect to the second derivatives 



Q" + A n Q + A 12 M + 7Q 3 
M" + A 21 Q + A 22 M 








(3.16) 



by the means of the transformation 



<1 
m 





-(3a 2 u 2 \f]5 



1 + f3a\ 1 + f3a\ 
The matrix A in (13.161) is symmetric A\ 2 = A 2 \ : 



Q 

M 



A 



ai + 



(X 2 UJ 2 yf]3 

1 + /3ai 



a 2 uj 2 \f]3 
1 + [3ai 1 + (3ai 



1 + /3ai 



(3.17) 



Instead of the second order system (I3.16P we will consider the following equivalent 4x4 
first order system 



(3.18) 



Obviously [0, 0, 0, 0] (the zero background) is the only equilibrium solution (the critical point) 
of the system. The pulse solutions are the trajectories of the system (13. 18j) that start and 
end at the equilibrium (homoclinic orbits). Thus, the investigation of solitary pulses is 
mathematically equivalent to studying homoclinic solutions. 
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4 Structure of solitary waves 

To investigate the structure of homoclinic solutions, we linearize the system in (13.181) near 
the critical point Q = M = Qi = Mi=0: 
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(4.19) 



The characteristic equation of the matrix A on the right-hand side is given by 

p 4 + (A u + A 22 ) p 2 + A n A 22 - A 21 A 12 = 
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Figure 1: The left figure shows the E = cross-section of the potential energy landscape 
U(Q,M) = 0. The Newtonian particle trajectory corresponds to a one-hump solution pre- 
sented in the right figure. 



Therefore, the values of p 2 coincide with the eigenvalues of the matrix 
that 

1 + Pen 

Thus, the condition 

uj\ — ai < 

makes det A < 0, causing A to have eigenvalues of opposite signs, which is a necessary 
condition for the existence of homoclinic orbits. Indeed, in the case of A having eigenvalues 
of the opposite signs, the 4x4 matrix A has two pure imaginary eigenvalues (square roots of 
the negative eigenvalue of —A) and one negative, and one positive eigenvalues. Therefore the 
nonlinear system has one- dimensional stable and unstable manifolds, and a two-dimensional 
center manifold (corresponding to the imaginary eigenvalues). 

It was first noticed in [l4| that the nonlinear system in (13.181) has Hamiltonian structure. 
If the kinetic, E, and potential, U, energies and the Hamiltonian, H, are introduce as follows 

E = i (Ql + Mf) , (4.22) 

U = A 12 QM + X - (A n Q 2 + A 22 M 2 ) + ^Q 4 , (4.23) 
H = E + U (4.24) 

then the system ( 13.181) takes the form 

d^Q l = -dH/dQ, d c M x = -dH/dM, 
d c Q = dH/dQx , d ( M = 8H/dM l . 



—A. It is easy to see 
(4.20) 

(4.21) 
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Soliton Q-M planeV=0.9 Soliton V=0.9 




Figure 2: The Newtonian particle trajectory (left) corresponds to the two-hump solitary 
wave (right). 



Since the Hamiltonian is a conserved quantity, d^H = 0, any trajectory issued from the 
critical point [0, 0, 0, 0] stays on the zero energy level surface H = for all time. Note that 
the surface H = is a 3D manifold in R 4 . The intersection of this 3D hypersurface with the 
hyperp lanes Qi = and Mi = is a curve T in the QM-plane 

U(Q, M) = A 12 QM + i (A n Q 2 + A 22 M 2 ) + ^Q 4 = (4.25) 

(the figure-eight shaped curve on the left in Fig. [1] and [2]). If for a given V there exits 
a homoclinic trajectory of (13.181) . then on this trajectory E + U = and since E > 0, 
necessarily U < 0. At the extrema of U its gradient is zero: 

8U BU 

— = A 22 M + A l2 Q = 0, — = A 12 M + A n Q + 7 Q 3 = 
dm oQ 

By eliminating M from the equations above, we obtain the cubic equation 

A 2 

l2 Q + A n Q + 7Q 3 = Q 



A 



22 



whose roots are easily found: 



-detA 



Q = 0, Q = ± x /^— = ± x /a 1 - uo 2 . 



Thus, VZ7 = at the points 



(0,0), (±^ ai -ujl ±^//?( ai -u;?)) 
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which are real if 

V 2 

thus producing the figure eight level curves. We already encountered this inequality above, 
see (14.2ip . After some algebra it can be rewritten as the following constraints on the traveling 
wave velocity: 

V <V<1, V = ^ujf/il + uf) (4.26) 
Thus, a possible velocity of the propagating pulse is bounded below. 



5 Numerical study of solitary waves 



The nonlinear system (I3.18P has the time reversal symmetry; therefore, if [Q, M, Qi, Mi](t) := 
u(t) is a homoclinic orbit, [Q, M, — Qi, — Mx](— t) also is (recall that Q\ and M\ are time 
derivatives of Q and M). 
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Distribuiton of solitons over the velocity spectrum 
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(a) Number of solitons per bin: 100 bins is de- (b) The hump distribution: number of humps vs. the 
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velocity. 

Figure 3: Statistics of solitary wave solutions. 



A priori it is not clear why any homoclinic solution would possess this symmetry, and it 
is quite likely that there exist non-symmetric homoclinic orbits; we plan to investigate them 
elsewhere. The characteristic property of a time reversal orbit is that at the symmetry point 
Qi = M\ = 0, and consequently the kinetic energy E must be zero; i.e., the symmetry point 
lies on the curve T, see (14.251) . Moreover at the symmetry point the trajectory is orthogonal 
to T (for an illustration, see the QM diagrams on the left of Fig. [T]and Fig. [2]). 

Our search algorithm for finding solitary-wave solutions is based on the following mini- 
mization idea. If for a given value of the propagation velocity V there exists a homoclinic 
orbit with the time-reversal symmetry, then at some point both the kinetic and potential 
energies are zero. The algorithm takes the initial condition u from a domain S on the 
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Figure 4: Solitary waves examples V = 0.84322 and V = 0.91461. The left and right figures 
illustrate four-hump and eight-hump solitary wave solutions. 



Energy Norm of Solitons Over the Velocity Spectrum 



zero-energy surface, near the critical point (0, 0, 0, 0) and in the direction close to that of the 
unstable eigenvector of the linearized problem. Then the following optimization problem is 
posed: Determine 

$(V) = min min E[u(C|u )] (5.27) 

unGS Co<C<Co+to 

where E is the kinetic energy, u(£|uo) is 
the solution of (13.181) with the initial con- 
dition u(0|uo) = uo, recall that u = 
[Q, M, Qi, Mi}. The parameter r is the ex- 
pected "width" of the pulse. Since E = 
at C = we take C > Co to obtain a non- 
trivial solution for the energy minimization 
problem. Computation of any particular 
value of min£ 0< £<£ 0+To -E[u(C|u )] involves a 
numerical solution of the nonlinear system 
of ODEs. 

The search of the optimal initial datum 
is stochastic and is organized via a version 



10* 



3JJ. On each step 



10' 




of simulated annealing 
the initial datum is obtained by sampling a 
random distribution with the density deter- 
mined by the results of the previous step (see 
3 Q for more detail). If $(V) = then 
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Figure 5: The energy per hump vs. the veloc- 
ity 



there exists a homoclinic solutions with velocity V. 
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Figure 6: Stable propagation of a eight-hump solitary wave; V = 0.822. 



When the kinetic energy possesses several local minima along the trajectory the corre- 
sponding solitary wave has the multi-hump structure. Figures [1] and [2] illustrate this phe- 
nomenon. The figure-eight shaped curves on the left correspond to the E = cross-section 
of the potential energy landscapes; the curves inside the domains represent the Newtonian 
particle trajectories in the QM configuration space. The graphs on the right show the 
profiles of the corresponding solitary wave solutions. Fig. [1] illustrates a typical one-hump 
solution. In contrast, the trajectory shown in Fig. [2] has a point of the nearest approach to 
the boundary where the kinetic energy attains a local minimum. The resulting solution has 
a two-hump structure. Multi-hump solutions correspond to more complicated trajectories. 
Each of these trajectories has the return point at which it has the normal incidence with the 
E = contour. 




Figure 7: The initial solitary wave pulse with a perturbation added (left); Evolution of this 
pulse governed by the PDEs (right). 



For the fixed set of physical parameter values, the shape of the potential energy landscape 
is controlled by the pulse velocity V via the coefficients a± and a<i in (13.151) . We investigated 
numerically the set QJ of values of V which give rise to homoclinic orbits; in some sense 
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one might think of these Vs as the "spectrum" of the problem. For numerous applications 
with soliton-like solutions the velocity value is known to change continually (a continuous 
spectrum). However, for the Maxwell-Duffing model under consideration our numerical 
investigation demonstrates that the spectrum 93 contains both an interval of a continuous 
spectrum and a discrete subset of parameter values V for which a wave solution exists. One 
of the principal issues is to understand the correspondence between types of solitary wave 
solutions and values of V G 93. 

We first investigated numerically the distribution density of the values of V, which give 
rise to homoclinic orbits. For all numerical computations of this section we adopted the 
following values of the nondimensional physical parameters: 



1, UJ2 = 5, 7 = 0.01, = 0.5 



(5.2* 



For uji — 1 the allowable range of values of V from ( 14.26ft is given by l/y/2 < V < 1. The 
plot in Fig. 3(a) illustrates the density distribution of V G 93 on the interval [0.73,0.95]. 
The search algorithm tested potential values of V on the grid 8V = 10~ 4 . The plot depicts 
the number of "successful" homoclinic orbits per velocity interval AV (a "bin"); in this 
particular case the value has been chosen as AV = 0.002. 




200 



180 



Figure 8: Solitary wave collisions: two-hump solitons with V = 0.9 and V = —0.9 (left); an 
eight-hump soliton with V = 0.89 and a phase- inverted soliton with V = —0.75 (right). 



Our numerical computations show that on a rather small interval 93 c = [0.73, 0.7642] at 
the low end of the spectrum, every attempt of computing a homoclinic orbit was successful 
(20 orbits per bin). These results stay consistent with the refinement of the computational 
grid size SV. 

All the solitary wave solutions in 93 c are of the one-hump variety; note, however that 
the single-hump solitons are not exclusively confined to the lower end of 93. Elsewhere the 
spectrum density is very low, and the solitons are mostly of a multi-hump kind. Somewhat 
arbitrarily, we define a hump as a local maximum of the electric field e, which is at least 
50% of the global maximum. 

Next we studied the distribution of the different type of solitary wave solutions on the 



interval of velocities [0.73, 0.95]. The figure (Fig. 3(b) ) gives a very clear idea of the placement 
of solitons according to the number of humps, which ranges from one to ten. Some typical 



11 



soliton profiles for four- and eight-hump solutions with V = 0.84322 and V = 0.91416 
respectively are collected in Fig. HI 




Figure 9: a) Initial pulses for propagation study; b) Energy dissipation for the given initial 
pulses. Larger Gaussian waves quickly shed energy while breaking up into near-solitary 
waves. The near solitary-waves slowly lose energy while converging to solitary waves. 

Different types of solitary wave solutions have different energy values. Because of the 
multi-hump nature of these solutions it is convenient to introduce the energy of the electro- 
magnetic field per one hump. We analyzed dependence of the electromagnetic field energy 
S per one hump versus velocity of solitary wave, see Fig. [51 Here S is defined as 



1 f°° 

S= 2Nj t^M+^ft*)] 



dt, 



where N is the number of humps. As follows from this figure, in the log-log coordinates the 
energy increase is very well approximated by a linear function. The least square fit of the 
data from this figure shows that the energy increases approximately as a polynomial of fifth 
degree in V. 

6 Formation, stability, and interaction of solitary waves: 
computer simulations 

In this section through direct numerical simulations we study evolution of waves as well as 
wave interactions. We consider formation of solitary wave solutions from arbitrary initial- 
boundary condition, stability of traveling waves under small perturbations and stability 
under strong perturbations due to wave collisions. In all numerical simulations of this section 
we use the same values of physical parameters (I5.28|) as in Sec. [51 
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Numerically we solve the signaling problem for (12.71) . In other words we give boundary 
conditions on either one or both ends of the spatial interval (0, L); as initial conditions we 
assign zero values for all the variables, which corresponds to propagation in a quiescent 
medium. For solving the initial-boundary value problem for the system in (12. 7j) we devised 
a simple fractional step numerical method. 

Because the first two equations in (12.71) are hyperbolic PDEs while the rest are ODEs 
the choice of the fractional steps is extremely natural: on the first half-step we propagate 
the PDE part of the governing equations, and on the second half-step we march according 
to the system of ODEs. The resulting ODE system is solved by using the midpoint rule, 
while the PDEs are solved by the explicit McCormack method 33[. The midpoint rule and 
the McCormack method are both second order accurate. To increase the accuracy of the 
fractional step method we utilize the Strang split approximation 34|, which results in the 
second order convergence of the final numerical scheme. 

For many of the solitary-wave solutions discussed in Sec. we ran direct numerical 
simulations on the model with these solitary waves as input pulses. All the waves tested, 
even the ones of a rather intricate shape, propagate with constant speed and without any 
shape distortion. See, for example, Fig. [6] where propagation of an eight-hump soliton is 
depicted. 




Figure 10: Evolution of the 15 exp(— O.lr 2 ) Gaussian. A "sharp" Gaussian quickly evolves 
into a near-solitary wave, leaving some disturbance in the wake. The speed of the near- 
solitary wave is significantly higher than the speed of propagation of the radiation; thus the 
wave quickly leaves the disturbance behind 

This suggests that the solitary waves are (nonlinearly) stable with respect to numerical 
perturbations. We remark that although because of the scale of Fig. El the pulses appear 
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rather singular, they are in fact completely smooth and numerically resolved. The numerical 
resolution of this computation is 8 mesh intervals per unit length, which provides about 60 
mesh points per each hump of the traveling wave. Similarly, fine computational meshes are 
employed in all the simulations below. 

The issue of stability can be addressed analytically by studying the linearization of the 
system of partial differential equations (12.71) about arbitrary traveling wave solutions and 
analyzing the corresponding linear evolution operator. Our analysis showed that this op- 
erator is skew-Hermitian in L 2 with the appropriate norm. Therefore the spectrum of the 
evolution o per ator is pure imaginary and the traveling wave solutions are neutrally linearly 
stable (see [28j for detail). 



Gaussian 15*exp(-.05*t ) 




Figure 11: Evolution of the 15 exp(— 0.05r 2 ) Gaussian. A medium size Gaussian evolves into 
two waves. The velocity of the smaller wave is on the order of the velocity of radiation. 

To further elucidate the issue of stability we consider stability with respect to a finite- 
size perturbation in the initial wave shape. This situation is illustrated in Fig. [7J To the 
two-hump numerical soliton we add a rather substantial perturbation and employ the thus 
obtained functions as boundary data for the system of partial differential equations (|2.7p . 
As the result of evolution, the solution relaxes to the solitary wave shape followed by a 
low-amplitude "continuous radiation" . 

Stability with respect to strong perturbations due to collision of two traveling wave 
solutions is illustrated in Fig. [HJ We take two solitary waves obtained by numerical solutions 
of ODEs and use these solutions as the boundary conditions for the PDEs. The left part 
of the Fig. [HJ shows collision of two-hump solitary wave solutions. The right part of the 
figure shows collision of eight-hump and one-hump solitary waves. In both cases collision of 
solitary waves leads to formation of the steady state solutions. The collisions are followed 
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Figure 12: Evolution of the 15 exp(— O.Olr 2 ) Gaussian. Large Gaussian quickly breaks up 
into four near-solitary waves, leaving some disturbance in the wake. The waves become 
more separated over the time, since the near solitary waves with higher amplitude have 
higher velocities. 



by emission of a small amplitude continuous radiation and a residual phase shift. 

A soliton nature of solutions of (12.71) is further confirmed by the set of numerical simu- 



lations we present next. We consider propagation of solutions with the pulses in Fig. 9(a) 



given as a series of boundary conditions at the x = boundary. The soliton of velocity 
V = 0.75 (see Fig. [I]) propagates in a stable fashion, while its least-squares approximation 
by a Gaussian 9.65 exp — 0. It 2 approaches the soliton shape after shedding a small amount 
of residual continuous radiation. These time evolutions are not included for space saving 
(the energy dissipation curves for these cases show conservation of electromagnetic energy, 
see Fig, gfc)] ). 

In the next three figures (Fig. [T0l- fT2|) we present evolutions of the larger Gaussian pulses 
from Fig. 9(a) Evolution of the sharpest Gaussian (a = Vo) is displayed in Fig. [10] Very 



fast the solution forms a solitary wave that moves with constant velocity with no shape 
change. It is followed by low magnitude oscillations whose leading edge also moves with 
constant speed. During the evolution, the oscillatory part disperses more and more. This 
part of the solutions appears to be of a nonlinear nature; it will be studied separately. We 
note that although because of the scale of the figure, the pulse appears very sharp, it is in 
fact completely smooth with "width" about 20 and about 150 computational mesh points 
within the pulse. 

The evolution of a wider Gaussian, a = V^IO (see Fig. [TTI) is similar with a very interesting 
distinction. Now the leading soliton is trailed by a slower low amplitude soliton. The latter 
is followed by low amplitude oscillations that again lag behind and disperse. The waves 
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become more separated over time because the solitary waves with higher amplitude have 
higher velocities. Finally, the widest Gaussian, a = 5^/2, develops into a train of four 
solitons, see Fig. [T2] 

To characterize the energy exchange between the propagating pulse and the medium, in 
Fig. 9(b) we present plots of the total electromagnetic energy as a function of time for all 
the input profiles from Fig. 9(a) For the soliton solution there is a dynamic equilibrium 
between the energy stored in the medium and the electromagnetic energy of the pulse. In 
case of the input impulse being not a soliton, the balance between the medium and the pulse 
is violated, which leads to the dissipation of the electromagnetic energy into the medium. 



7 Concluding Remarks 

In this paper we considered propagation of extremely short pulses in a nonlinear medium, 
which is characterized by both electric and magnetic resonance responses. Interaction of 
the electromagnetic field with the medium was described in the framework of the Maxwell- 
Duffing model. In particular we employed the classical Maxwell-Lorenz model for describing 



the magnetic resonance, [15] . For describing the interaction of the electric field component 
with the medium we used a generalized Maxwell-Lorentz model which takes into account 
cubic anharmonism of the polarization response (i.e., the Maxwell-Duffing system). Our 
findings demonstrate that the model supports a wide array of traveling-wave solutions. We 
investigated the structure and properties of these solutions through a combination of anal- 
ysis and numerical modeling. We determined that the family of traveling-wave solutions is 
parameterized by one parameter, which is the velocity of a steady wave solution, normalized 
by the speed of light in vacuum. The spectrum 53 contains both an interval of a continu- 
ous spectrum and a discrete subset of parameter values for which a traveling-wave solution 
exists. Computer modeling demonstrated a multi-hump structure of these solutions. Their 
multi-hump nature suggests to characterize solitary wave solutions by a number of humps 
(types). All types are determined by not overlapping sets of velocities. 

Direct numerical simulations showed that solitary-wave solutions are dynamically stable. 
This dynamical stability is consistent with the analysis of the system linearized about solitary 



wave solutions 28]. Stability of these solutions with respect to strong perturbations was 
studied by means of solitary wave collisions. Computer simulations indicated nearly elastic 
nature of scattering followed by a residual excessive radiation and a phase shift. In addition 
to traveling-wave solutions, numerical simulations demonstrated presence of another type of 
nonlinear oscillatory solutions with extended tail. 
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